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ABSTRACT 


Several  techniques  for  estimation  of  the  power  spectral  density  (PSD),  or 
simply  the  spectrum,  of  discretely  sampled  short  data  sequences  are 
investigated  in  this  thesis.  These  techniques  are  directly  based  on  the  fast 
Fourier  transform  (FFT)  of  the  data  sequence.  However,  an  extrapolated 
version  of  the  original  data  sequence  is  used  instead  of  the  original.  The  aim 
is  to  enhance  the  resolution  of  the  true  spectral  components  of  the  signal 
without  adding  any  false  frequency  components.  Several  data  extrapolation 
models  are  analyzed  and  their  performance  at  different  signal  to  noise  ratios 
(SNRs)  and  model  orders  are  examined. 
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I.  INTRODUCTION 


This  thesis  investigates  several  techniques  for  enhancing  the  spectrum 
estimation  of  short  data  sequences.  A  short  data  sequence  is  defined  as  one  in 
which  the  spectral  resolution  required  is  of  the  same  order,  or  better,  as  the 
reciprocal  sequence  length.  The  techniques  to  be  investigated  will  involve 
extrapolating  the  original  data  sequence  and  then  performing  a  basic 
periodogram  on  this  new  data  sequence.  By  extending  the  original  data 
sequence  we  hope  to  enhance  the  true  spectral  components  and  the  spectral 
resolution  of  the  signal  without  adding  any  false  spectral  components. 

The  method  of  spectrum  estimation  to  be  used  on  all  extended  data 
sequences  is  the  basic  periodogram.  The  periodogram  is  defined  as 


where  x[n]  is  the  extrapolated  data  sequence  and  M  is  the  number  of  points  of 
the  extended  data  sequence. 

This  thesis  is  arranged  into  5  chapters  and  4  appendices.  Chapter  II  is 
devoted  to  data  extrapolation  via  linear  predictive  filtering.  In  Chapter  HI  the 
data  is  modeled  as  the  impulse  response  of  an  infinite  impulse  response  (HR) 
filter.  Data  extrapolation  by  dominant  eigenvectors  is  discussed  in  Chapter 
IV.  In  the  last  chapter  some  general  conclusions  of  the  work  carried  out  in 
this  thesis  are  presented,  and  suggestions  for  future  investigations  are  given. 
Appendix  A  covers  several  other  methods  of  AR  parameter  estimation. 
Appendix  B  gives  some  results  from  a  causal  Prony's  method  of 
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extrapolation.  Appendix  C  contains  some  results  from  some  miscellaneous 
test  sequences  while  computer  simulation  programs  are  included  in 
Appendix  D. 
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II.  LINEAR  PREDICTIVE  HLTERING 


A.  DEHNUION 

In  this  chapter,  all  original  test  data  sequences  will  be  extrapolated  by 
linear  prediction.  In  linear  prediction,  a  given  set  of  sample  values  of  a 
random  sequence  x  is  used  to  predict  the  values  of  the  sequence  some  time 
into  the  futiue.  We  assume  that  only  P  previous  values  of  the  sequence  are 
used  in  the  prediction.  Now  the  problem  is  to  find  the  prediction  coefficients 
-a^,  -a2,  ...,  -ap  to  produce  an  estimate  of  the  current  value  of  the  random 
process  x[n]  from  P  previous  values.  The  estimate  is  written  as 

i[n]  =  -a\x[n  - 1]  -  -  2]-.  ..-d^n  -  P]  (Zl) 

where  *  represents  the  complex  conjugation. 

The  error  in  the  estimate  is  given  by  the  difference 

in]  =  x[n]  -  x[n]  =  (2-2) 

k=0 

where  ao  is  defined  as  1. 

The  approach  is  to  choose  the  prediction  filter  coefficients  to  minimize 
the  mean-square  error 

o?  =  E[le[nf  ]  =  E[|x[n]-x[nf  ].  (Z3) 

To  find  the  optimal  prediction  filter  coefficients  we  apply  the  Orthogonality 
Theorem  [Ref.  l:p.  308  ]  which  states  that 

E[x[n-i]e*[n]]  =  0  for  i  =  1,2,...P  (Z4) 


or 
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p 

=  fori  =  l,2,...,P  (2.5) 

/=l 

where  fxx  is  the  autocorrelation  function  of  x. 

The  minimum  prediction  error  is  found  to  be 

P 

k=\ 

These  equations  are  identical  to  the  Yule-Walker  equations  for  autoregressive 
(AR)  processes  [Ref.  2;p.  157].  Therefore,  the  optimal  linear  prediction 
coefficients  are  just  the  AR  parameters,  when  the  order  of  the  AR  process  is 
identical  to  the  order  of  the  linear  predictor. 

B.  AR  PARAMETER  ESTIMATION 

For  good  estimates  of  the  AR  parameters  or  prediction  coefficients,  a 
maximum  likelihood  estimator  (MLE)  is  usually  employed.  Several  methods 
were  tried  in  this  thesis  including  the  autocorrelation  method,  the  covariance 
method,  the  forward-backward  method  (modified  covariance  method),  the 
modified  forward-backward  method,  and  the  Burg  method.  The  two 
methods  which  proved  to  be  the  most  robust  were  the  forward-backward 
method  and  the  modified  forward-backward  method.  Explanations  and 
results  from  the  other  methods  are  given  in  Appendix  A. 

1.  Forward-backward  Method  (Modified  Covariance) 

For  an  AR  process  the  optimal  forward  predictor  is  given  by  Eq.  (2.1). 
The  corresponding  optimal  backward  predictor  is  given  by  [Ref.  2:p.  226] 

P 

x[n]='^-af^x[n+k].  (Z7) 

jk=l 

The  forward-backward  method  estimates  the  prediction  coefficients  by 
minimizing  the  average  of  the  estimated  forward  and  backward  prediction 
error  powers.  The  average  of  the  estimated  forward  and  backward  prediction 
error  powers  is  given  by 


4 


1 


p  =  i(p/+p'’)  (Z8) 

where  the  forward  and  backward  prediction  error  powers  are  respectively 
defined  as 


,  N-l 

y 

NT-P  " 
^  ^  n=P 


4”1+  '^a[^]x[n~k\ 
k=\ 


N-P 


N-l-P 

I 

n=0 


x[n]+  ^a*[k]x[n  +  k] 
k=l 


(2.9) 


and  N  is  the  number  of  data  points  in  the  random  sequence  x.  To  minimize 
Eq.  (2.8)  we  differentiate  p  with  respect  to  the  real  and  imaginary  parts  of  a[k] 
for  k=l,2,...P  and  set  the  expression  equal  to  zero.  After  simplification  this 
becomes 


P  /N-l  N-l-P  ^ 

'^x[n~k]x*[n-i]+  +  + 

k=\  \n~P  n=0 

/N-l  N-l-P  ^ 

=  -  '^x*[n]x[n  +  t] 

<n=P  n=0  > 


(2.10) 


for  /=1,2,...,P. 

Defining 


/N-l 


N-l-P 


,  (2-11) 


Eq.  (2.10)  can  be  written  in  matrix  form  as 


■c„[l,l]  C„[l,2]  -  C^[1,P)' 

■c„[i,o]' 

Cxx[2.1]  C„[2,2]  -  C„[2,P] 

42] 

=  — 

^jtx[2.0] 

C:„[P,1]  C„[P,2]  -  C„[P,P), 

API 

.CxxlP.O]. 

(212) 
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The  first  matrix  is  a  P-by-P  autocorrelation  matrix  while  the  last 
matrix  is  called  the  cross  correlation  vector.  The  prediction  coefficients  are 
easily  solved  for  from  Eq.  (2.12). 

Intuitively,  the  order  P  should  be  as  large  as  possible  in  order  to  have 
a  large  data  base  for  the  predictor.  The  use  of  too  large  a  value  of  P  can  give 
rise  to  spurious  spectral  peaks  in  the  AR  spectral  estimation  using  the 
modified  covariance  method  (Ref.  3;p.  347].  For  the  best  performance  of  the 
forward-backward  linear  predictor  (FBLP)  method,  Lang  [Ref.  4:p.  720]  suggests 
the  value 


where  once  again  N  is  the  original  data  length. 

Akaike  [Ref.  5:pp.  203-217]  proposes  two  other  methods  of  model 
order  selection  based  on  the  estimated  prediction  error  power. 

a.  Test  Data 

Two  types  of  test  data  were  used  for  computer  simulation  of  all 
methods.  The  first  set  was  obtained  from  a  paper  by  S.  M.  Kay  and  S.  L. 
Marple  [Ref.  6:p.  1380-1414].  The  data  is  made  up  of  three  sinusoids  and  a 
colored  noise  process  which  is  obtained  by  filtering  a  white  Gaussian  process. 
The  three  sinusoids  are  at  fractional  frequencies  of  the  sampling  frequency  of 
0.10,  0.20,  and  0.21  and  have  SNR's  of  +10,  +30,  and  +30  dB  respectively.  SNR 
is  defined  as  the  ratio  of  the  sinusoid  power  to  the  total  power  in  the 
passband  noise  process.  The  noise  process  passband  is  centered  at  0.35  and  has 
a  bandwidth  of  0.30.  Figure  1  is  a  plot  of  this  64  point  test  sequence  and  the 
true  PSD.  Figure  2  is  the  periodogram  of  the  original,  nonextrapolated  data 
sequence  using  a  Hamming  data  window.  A  Hamming  window  is  used  in  all 
of  the  periodogram  estimations  in  this  thesis. 
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The  second  type  of  test  data  is  made  up  of  two  sinusoids 
corrupted  by  white  noise  .  This  is  also  a  64  point  test  sequence.  This  type  of 
test  data  will  be  composed  of  different  frequencies,  phases,  and  SNRs  to  see 
how  each  method  of  spectrum  estimation  resolves  the  sinusoidal  frequencies. 
A  sampling  frequency  of  1000  Hz  is  used  to  generate  the  signals,  therefore  the 
maximum  frequency  of  any  single  sinusoid  is  500  Hz.  Using  ordinary  FFT 
routines  on  non-extended  64  point  data  sequences,  the  spectral  resolution  can 
be  as  good  as  0.8(1 /64msec)=l  2.5  Hz.  Using  extended  data  sequences,  it  will  be 
shown  that  the  resolution  can  be  improved  to  2.0  Hz  for  an  appropriate  SNR. 

For  comparison  purposes,  the  periodograms  of  the  extrapolated 
data  sequences  are  plotted  with  the  autoregressive  (AR)  spectrum  estimations 
of  the  same  data  sequences.  The  AR  spectrum  is  defined  as 


(2.13) 


where 


(2.14) 


and  of  course  the  AR  coefficients  are  identical  to  the  prediction  filter  coefficients. 

The  dashed  lines  on  each  plot  represent  the  positions  of  the  true 
frequencies  of  the  test  data. 
b.  FBLP  Results 

Figure  3  is  the  extrapolated  version  of  the  Kay  and  Marple  (KM) 
data.  The  data  was  extrapolated  by  using  linear  prediction  as  in  Eq.  (2.1).  The 
prediction  coefficients  were  obtained  from  the  FBLP  method  utilizing  Eq. 
(2.12).  The  data  length  was  extrapolated  to  320  points  using  a  model  order  of 
16. 
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Care  must  be  taken  when  extrapolating  the  data  sequence.  If  the 
correlation  matrix  in  the  calculation  of  the  prediction  coefficients  is  not 
positive  definite,  the  resulting  poles  in  the  filter  of  the  AR  model  are  not 
guaranteed  to  lie  within  the  unit  circle.  This  could  cause  the  extrapolation  to 
grow  without  bound  and  will  have  devastating  effects  on  the  spectrum 
estimation.  The  forward-backward  method  does  not  guarantee  a  stable  all 
pole  filter,  although  in  practice  the  extrapolation  is  usually  boimded.  We 
should  check  to  ensure  all  poles  of  the  AR  process  lie  within  the  unit  circle 
before  we  extrapolate  the  data. 

Figure  4  shows  the  periodogram  of  the  extrapolated  KM  data  and 
the  respective  AR  spectral  estimation.  Both  methods  correctly  calculate  the 
three  sinusoidal  frequencies  present.  The  colored  noise  process  is  also  given  a 
fair  representaion.  When  the  data  is  extrapolated  to  576  points,  an  even  better 
PSD  estimate  is  obtained  from  the  FBLP  extrapolation.  Figure  5  shows  the 
results. 

For  the  two  sinusoids  corrupted  with  white  noise,  good  results 
are  obtained.  Figure  6  shows  the  spectrums  of  two  sinusoids  with  a  10  Hz 
separation  and  a  SNR  of  12  dB.  The  original  data  sequence  was  extrapolated 
to  320  points.  Both  sine  waves  have  a  0  degree  initial  phase  angle.  The 
frequencies  are  correctly  calculated  using  the  periodogram  of  the  FBLP 
extrapolated  sequence  while  the  AR  spectral  estimation  is  unable  to 
distinguish  the  two  closely  spaced  sinusoids. 

When  the  initial  phase  angle  between  the  two  sine  waves  is 
changed  to  45  degrees,  the  periodogram  of  the  FBLP  extrapolated  sequence. 
Figure  7,  distinguishes  the  two  frequencies,  even  though  one  of  the  estimated 
spectral  components  is  slightly  off.  The  AR  spectrum  estimation  is  unable  to 
resolve  them.  Once  again  the  data  was  extrapolated  to  320  points. 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


Figure  4.  Periodogram  of  Extrapolated  KM  Data  (320  points)  and  Respective 

AR  Spectrum  Estimation 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


Figure  7.  Periodogram  of  Extrapolated  Sequence  and  Respective  AR 
Spectrum  Estimation  of  two  Sine  Waves  (45°  Initial  Phase  Angle) 
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When  the  initial  phase  angle  between  the  two  sine  waves  is 
changed  to  90  degrees,  in  essence  a  signal  consisting  of  a  sine  wave  and  a 
cosine  wave  with  frequencies  that  are  10  Hz  apart,  neither  method  is  able  to 
distinguish  the  two  frequencies. 

2.  Modified  Forward-Backward  Method 

The  modified  forward-backward  method  is  a  two-step  procedure  for 
modifying  the  FBLP  method  [Ref.  3:p.  347].  First,  we  compute  the  P 
eigenvalues  of  the  correlation  matrix  given  in  Eq.  (2.12).  To  do  this,  we  use 
the  singular  value  decomposition  (SVD)  which  is  applied  to  the  2(N-P)-by-P 
data  matrix  A  given  by 


■  x[P] 
x[P--l] 


x[P  +  l] 
x[P] 


x[N-l]  x*[2]  x*[3] 
x[N-2]  x*[3]  x*[4] 


z[l]  jc[2]  •••  x[N-P]  X •[?  +  !]  x*[P  +  2] 


x*[N-P-i-l]‘ 
:c*[N-P  +  2] 

x*[N] 

(2.15) 


where  N  is  the  length  of  the  original  data  sequence.  Note  that  the  correlation 
matrix  in  Eq.  (2.12)  is  equal  to  A^^A. 

Once  the  P  eigenvalues  are  computed,  we  divide  the  data  space  into 
two  parts:  (i)  the  signal  subspace  spanned  by  the  eigenvectors  associated  with 
the  K  largest  eigenvalues,  and  (ii)  the  noise  subspace  spanned  by  the 
remaining  (P-K)  eigenvectors.  The  notion  used  is  that  the  K  principal 
eigenvectors  of  the  correlation  matrix  are  due  to  the  signal  components  amd 
are  less  sensitive  to  perturbations  caused  by  noise  than  the  remaining 
eigenvectors.  In  the  case  of  noiseless  data  consisting  of  K  sinusoids,  the 
correlation  matrix  has  K  nonzero  eigenvalues  and  (P-K)  zero  eigenvalues. 
Thus  in  a  practical  case  of  noisy  data,  by  keeping  the  K  largest  eigenvalues  of 
the  correlation  matrix  and  ignoring  the  rest,  we  are  in  effect  attempting  to 
increase  the  SNR  by  approaching  the  noiseless  case. 
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Let  v  1,1)2,  •••,  i)K  denote  the  K  largest  eigenvalues  and  ci,  62, ck 
denote  their  associated  eigenvectors,  respectively.  We  use  these  eigenvalues 
and  eigenvectors  to  compute  the  following  P-by-1  estimate  of  the  prediction 
coefficients: 


where 

b  =  [x[P  +  llx[P  +  2],...,jc[N],x  •[l],x  *[2] X  *  [N  -  P]].  (2.17) 

For  the  modified  FBLP  (MFBLP)  method,  we  have  removed  the 
imdesirable  effects  of  the  noise  subspace  eigenvectors,  therefore  we  may 
increase  the  predictor  order  P.  Tufts  and  Kumaresan  [Ref.  7:pp.  975-989]  have 
experimentally  determined  the  value  of  the  predictor  order  to  be 

P  =  3N/4  (117) 

where  N  is  the  original,  nonextrapolated  data  length. 

Since  the  computer  simulation  uses  64  data  points,  a  model  order  of 
48  is  selected. 

a.  Modified  FBLP  Results 

Many  test  runs  are  made  on  the  signal  consisting  of  2  sine  waves 
corrupted  by  white  noise.  Sinusoidal  frequencies,  phases,  and  SNRs  are 
changed  for  comparison  purposes.  Table  1  specifies  the  parameters  for  each 
run  and  the  corresponding  figure  number.  In  each  test  case,  the  data  sequence 
is  extrapolated  to  a  data  length  of  320  points  using  linear  prediction  with  the 
prediction  coefficients  obtained  from  the  MFBLP  method.  Figures  8  through 
24  show  the  results. 


TABLE  1.  PARAMETERS  FOR  TEST  SEQUENCES 


CASE# 

SINUSOIDAL 

FREQUENCIES 

(HZ) 

SNR  (DB) 

RELATIVE 
INITIAL  PHASE 

FIGURE  # 

1 

100  and  107 

12 

0° 

8 

2 

100  and  107 

9 

0® 

9 

3 

100  and  107 

6 

0° 

10 

4 

100  and  107 

12 

45® 

11 

5 

100  and  107 

9 

45® 

12 

6 

100  and  107 

6 

45® 

13 

7 

100  and  107 

3 

45® 

14 

8 

100  and  107 

12 

90® 

15 

9 

100  and  107 

9 

90® 

16 

10 

100  and  107 

6 

90® 

17 

11 

100  and  107 

3 

O 

o 

18 

12 

100  and  105 

12 

0® 

19 

13 

100  and  105 

12 

45® 

20 

14 

100  and  105 

18 

45® 

21 

15 

100  and  105 

12 

90® 

22 

16 

100  and  105 

18 

90® 

23 

17 

100  and  102 

37 

0® 

24 

In  all  test  cases,  the  periodogram  of  the  MFBLP  extrapolated 
sequence  performs  as  well  or  better  than  the  AR  spectrum  estimation.  Figure 
8  shows  that  even  though  both  methods  correctly  estimate  the  frequencies, 
the  periodogram  of  the  MFBLP  extrapolated  sequence  sometimes  has  better 
defined  spectral  peaks.  In  case  10,  depicted  in  Figure  17,  the  periodogram  of 
the  MFBLP  extrapolated  sequence  correctly  distinguishes  the  two  frequencies 
while  the  AR  spectrum  estimation  is  unable  to  do  so.  Even  in  test  case  17 
where  the  two  sine  waves  are  only  2  Hz  apart,  the  periodogram  of  the  MFBLP 
extrapolated  sequence  has  well  defined  peaks,  as  shown  in  Figure  25,  while 
the  AR  spectrum  estimation  just  barely  distinguishes  the  two  frequencies. 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


MODIFIED  FBLP  EX  TENSION,  ORDER  48,  SNR=9,  PM=45 


AR  SPECTRUM  ESTIMATION 


Figure  12.  Case  5  Spectrum  Estimates 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


Figure  16.  Case  9  Spectrum  Estimates 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


Figure  19.  Case  12  Spectrum  Estimates 
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Figure  20.  Case  13  Spectrum  Estimates 
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Figure  23.  Case  16  Spectrum  Estimates 


34 


1 


When  the  MFBLP  extrapolation  technique  is  tested  on  the  KM 
sequence,  the  periodogram  of  the  extrapolated  sequence  is  considerably  better 
than  the  AR  spectral  based  estimation.  In  general,  the  modified  FBLP  method 
utilizes  prior  knowledge  of  the  number  of  sinusoids  present  in  the  signal. 
When  the  K  largest  eigenvalues  are  chosen,  we  are  assuming  that  there  are  K 
sinusoids  present  in  the  signal.  If  one  of  these  K  largest  eigenvalues  is 
associated  with  the  noise,  inaccurate  estimates  will  be  computed.  This  can  be 
seen  in  Figure  26.  Because  the  KM  data  is  made  up  of  3  sinusoids,  the  3 
largest  eigenvalues  and  corresponding  eigenvectors  are  chosen.  However  in 
this  case,  one  of  these  eigenvalues  is  associated  with  the  noise  so  the  signal  at 
the  0.1  fractional  frequency  is  not  detected  in  the  AR  spectral  estimation. 
However  when  the  same  prediction  coefficients  are  used  in  the  data 
extrapolation,  the  resulting  periodogram  is  a  more  accurate  spectral 
estimation  than  the  AR  spectral  estimation.  Obviously,  the  extrapolation 
approach  is  more  forgiving  if  the  exact  number  of  sinusoids  is  unknown.  In 
Figure  27  the  four  largest  eigenvalues  and  corresponding  eigenvectors  were 
chosen,  one  more  than  the  actucd  number  of  sinusoids  present.  Once  again 
the  AR  spectral  estimation  is  incorrect,  picking  up  more  noise,  while  the 
periodogram  of  the  modified  FBLP  extrapolation  gives  a  reasonable  spectral 
estimate.  The  KM  sequence  is  extrapolated  to  a  length  of  320  points  before  the 
periodograms  are  computed. 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


Figure  26.  KM  Spectnun  Estimates  Using  the  Three  Dominant  Eigenvectors 
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RELATIVE  PSD  (dB)  RELATIVE  PSD  (dB) 


Figure  27.  KM  Spectrum  Estimates  Using  the  Four  Dominant  Eigenvectors 
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III.  PRONYS  METHOD 


A.  DEHNinON 

The  second  method  of  sequence  extrapolation  involves  modeling  the 
data  as  the  impulse  response  of  a  linear  time  invariant  filter.  This  approach 
is  shown  in  Figure  28.  Note  that  B(z)  has  q  zeros  and  A(z)  has  p  poles. 


Figure  28.  Impulse  Response  Estimation 


We  want  to  minimize  the  error  e[n]  in  some  sense.  If  we  would 
minimize  leads  to  non-linear  equations  and  is  very 

difficult  to  solve.  Therefore  we  will  use  an  indirect  method  to  solve  this 
problem.  The  equation  for  the  filter  is 

z[n]  +  a\x{n  - 1]+. .  .+apx[n  -p]  =  bo5{n)  +  b\5{n  - 1)+. .  -  q).  (3.1) 

The  coefficients  ai  and  bi  are  unknown  in  Eq.  (3.1)  and  we  want  the  output 
estimate  x[n]  as  dose  to  x[n]  as  possible. 
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Assume  that  we  have  N  data  points  and  that  N  ^  p+q+1.  The  transfer 
function  H(z)  of  the  filter  is 

.-1 


Mz)  l  +  fli2“'+...flp2“P 


We  can  write  Eq.  (3.1)  in  matrix  form  as 


x[0] 

0 

0  —  o' 

2[1] 

x[Q] 

0  —  0 

1 

fll 

. 

02 

x[q^\] 

.  a^/y  +  l-p] 

■ 

4N-1] 

x[N-2] 

.  jc[N-l-p] 

^0 

h 

h. 


or  more  compactly 


Xb' 

b' 

a  — 

0 

(3.2) 


(3.3) 


(3.4) 


where  0  =  [0, 0, 0]. 

For  the  case  where  N=p+q+l,  we  can  write  the  lower  set  of  equations  as 

X^a  =  0  (3.5) 


and  solve  directly  for  a.  Once  we  know  a  we  can  solve  the  upper  equations 
directly  for  b  where 

Xga  =  b.  (3-6) 

This  method  is  known  as  Fade's  approximation  [Ref  l:p.  4951. 

For  the  more  general  case  where  N  >  p+q+1  we  still  have  Eqs.  (3.5)  and 
(3.6),  however  now  Eq.  (3.5)  is  an  overdetermined  set.  If  we  solve  Eq.  (3.5)  in  a 
least  squares  sense,  then  we  can  still  find  b  from  Eq.  (3.6).  Solution  of  this 
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least  squares  problem  [Ref.  l:pp.  474-495]  leads  to  the  least  squares  normal 
equations 

(x7x^)a  = 

where 

This  is  known  as  Prony's  method  for  time  domain  IIR  filter  design.  This 
method  was  used  to  model  the  test  data  as  the  impulse  response  of  an  IIR 
filter.  This  impulse  response  is  then  extrapolated  over  a  longer  duration  and 
then  its  periodogram  is  computed. 

B.  NONCAUSAL  EXTRAPOLATION 

Figure  29  shows  the  original  Kay  and  Marple  data  sequence  and  the  64 
point  Prony's  impulse  response  estimate  of  this  data  sequence.  The  order  of 
the  numerator  and  denominator  in  the  transfer  function  of  Eq.  (3.2)  was 
chosen  to  be  16.  Several  model  orders  were  tried  before  settling  on  the  order 
of  16.  If  too  low  an  order  is  chosen,  there  is  not  enough  resolution  in  the 
periodogram  of  the  extrapolated  sequence.  If  too  high  an  order  is  chosen, 
spurious  spectral  peaks  arise.  Examples  of  varying  model  orders  are  given  in 
Appendix  B.  As  can  be  seen  in  Figure  29,  Prony’s  impulse  response  estimate 
of  the  data  sequence  for  an  order  of  16  is  a  very  good  approximation. 

By  choosing  the  appropriate  impulse  duration,  the  impulse  response  can 
be  extrapolated  as  far  as  we  like.  The  impulse  response  can  also  be 
extrapolated  in  both  forward  and  backward  directions.  We  can  extend  it 
backwards  by  reversing  the  original  data  sequence  and  performing  the  same 
computations  as  one  normally  would  for  the  forward  case.  Once  this  new 
Prony's  estimate  has  been  computed,  we  reverse  the  extension  and  attach  it  to 
the  front  of  the  original  data  sequence.  This  gives  a  noncausal  extrapolation 
of  the  original  data  sequence. 
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Figure  29.  Original  KM  Data  and  Frony's  Impulse  Response  Estimate  of  KM 

Data 
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This  was  the  method  performed  on  all  test  sequences.  A  strictly  causal  or 
forward  extrapolation  was  tried  but  the  results  were  not  as  good  as  the 
noncausal  extrapolation.  Some  examples  of  the  causal  extrapolation  are 
shown  in  Appendix  B.  The  length  of  the  extrapolation  is  provided  on  each 
figure. 

C  RESULTS 

For  the  KM  data,  the  noncausal  extrapolation  method  gave  excellent 
results.  Figure  30  shows  the  extrapolated  time  series  and  the  resulting 
periodogram.  The  three  sinusoids  present  are  correctly  identified  and  the 
colored  noise  passband  process  is  properly  represented. 

Figure  31  shows  the  extrapolated  sequence  and  corresponding 
periodogram  of  a  test  sequence  consisting  of  two  sinusoids  embedded  in 
white  noise.  The  two  sinusoids  have  frequencies  which  are  10  Hz  apart 
(fs  =  1000  Hz)  and  both  have  zero  initial  phase.  The  SNR  is  15  dB.  The 
periodogram  of  the  extrapolated  sequence  correctly  identifies  the  two 
frequencies  present.  When  the  initial  phase  of  one  of  the  sine  waves  is 
changed  to  45  degrees,  the  periodogram  of  the  extrapolated  sequence  once 
again  identifies  two  frequencies,  however  one  of  the  frequency  estimates  is 
slightly  off.  Figure  32  shows  the  results.  Figure  33  shows  the  method  does 
not  distinguish  the  two  frequencies  when  one  of  the  initial  phase  angles  is 
changed  to  90  degrees. 

Figures  34  through  36  show  the  extrapolated  time  series  and  resulting 
periodograms  for  two  sine  waves  5  Hz  apart  in  frequency  with  different  initial 
phase  angles.  The  SNR  is  30  dB  in  all  three  cases.  The  periodograms  in  each 
of  the  three  cases  correctly  distinguished  the  two  frequencies  present  when 
the  initial  relative  phases  were  0  degrees  (Figure  34),  45  degrees  (Figure  35), 
and  90  degrees  (Figure  36). 
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magnitude  in  dB  amplitude 


Figure  30.  Extrapolated  KM  Data  and  Respective  Feriodogram 
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Phase  Angle^  SNR  =  15  dB,  and  Resulting  Periodogram 


47 


PRONYS  EXTRAPOLATED  ESTIMATE  -  ORDER  16  NON-CAUSAL 


magnitude  in  dB  amplitude 


Figure  34.  Extrapolated  Sequence  of  two  Sine  Waves,  5  Hz  Apart,  0°  Initial 
Phase  Angle,  SNR  =  30  dB,  and  Resulting  Periodogram 
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magnitude  in  dB  amplitude 


Figure  35.  Extrapolated  Sequence  of  two  Sine  Waves,  5  Hz  Apart,  45°  Initial 
Phase  Angle,  SNR  =  30  dB,  and  Resulting  Periodogram 
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magnitude  in  dB  amplitude 
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siinipling  frequency 


Figure  36.  Extrapolated  Sequence  of  two  Sine  Waves,  5  Hz  Apart,  90°  Initial 
Phase  Angle,  SNR  =  30  dB,  and  Restdting  Feriodogram 
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The  periodogram  resulting  from  Prony’s  method  of  data  extrapolation  is 
unable  to  distinguish  the  two  sinusoids  which  are  5  Hz  apart  in  a  SNR  less 
than  approximately  27  dB,  Prony’s  method  of  data  extrapolation  resulted  in 
periodograms  which  were  not  as  accurate  as  the  MFBLP  method  in  lower 
SNRs.  In  high  SNRs  (20  dB  or  higher)  good  results  are  obtained  but  the 
periodograms  resulting  from  the  data  extrapolation  using  the  MFBLP 
prediction  coefticients  still  provide  better  estimates. 
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IV.  EIGENVECTOR  EXTRAPOLATION 


A.  DEHNmON 

In  this  chapter  we  develop  a  series  expansion  for  a  random  vector  in 
terms  of  an  arbitrary  set  of  orthonormal  vectors  [Ref.  8:pp.  261-262].  Once  this 
set  of  orthonormal  vectors  (eigenvectors)  are  computed,  we  will  use  the 
dominant  ones  to  extrapolate  the  data  sequence  and,  hopefully,  improve  the 
spectral  resolution. 

First  consider  a  set  of  N  orthonormal  vectors  qi,  q2, ...,  qN  that  satisfy  the 
following  two  requirements: 


q^qi 


1,  k  =  i 

0,  k*i 


(4.1) 


Let  X  denote  an  N-by-1  random  vector  with  zero  mean  and  correlation  matrix 
Rx.  The  requirement  is  to  express  the  random  vector  x  as  a  linear 
combination  of  the  set  of  orthonormal  vectors  qi,  q2, ...,  qN  as  follows: 

N 

X  =  S®iq«-  (4*2) 

1=1 

Equations  (4.1)  and  (4.2)  define  the  discrete  time  version  of  the  Karhunen- 
Lo6ve  expansion  [Ref.  8:p.  261].  The  coefficients  of  the  expansion,  Ci,  are 
defined  by  the  inner  product  between  the  vectors  x  and  qi,  given  by 

a,  =q[x 

=  x^q,-  for  i  =  1,2,..., N.  (4.3) 

The  coefficients  ai  are  obviously  random  variables  and  we  would  like  to 
choose  the  qi  so  that  these  coefficients  are  mutually  uncorrelated. 
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Because  the  random  vector  x  has  zero  mean,  we  observe  that  all  the 


coefficients  ai  also  have  a  mean  of  zero,  that  is 

E[a,]  =  0  fori  =  l,2,...,N.  (4.4) 

Therefore,  for  the  coefficients  of  expansion  to  be  uncorrelated,  we  require  that 
they  satisfy  the  condition 

covfajtttj]  =  E[(ajt  -E[ajt])(a,  -£[a,])] 


=  E[«itai] 

Vi,  k  =  i 

0,  k^i 

where  the  Ui  are  constants  to  be  specified.  From  Eq.  4.3  we  can  rewrite  Eq.  4.5 
as 


cov[ajt«,]  =  E[«]tai] 

=  E[q[xx^q,] 

=  qrE[«^]qi  (4-6) 

Since  by  definition  the  correlation  matrix  of  the  random  vector  x  equals 

R,  =  £{xx’'],  (4.7) 


we  can  simplify  Eq.  4.6  as 

covfajta,]  =  qfRxq,. 

Therefore  from  Eq.  4.5  and  Eq.  4.8  we  have 


q[Rxqf  = 


[vi,  k  =  i, 

lO  k*i’ 


(4.8) 


(4.9) 


If  we  multiply  both  sides  of  Eq.  4.1  by  Ui,  we  have 
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T  ^  ~ 


(4.10) 


Therefore  for  Eq.  4.9  to  hold  for  all  choices  of  k  and  a  particular  i,  it  is 
necessary  that  the  matrix  product  equal  Uiqi,  as  shown  by 

=  mi  (411) 

where  the  scalar  Uj  is  the  eigenvalue  of  the  correlation  matrix  Rx  and  the 
vector  q,  is  the  eigenvector  associated  with  the  eigenvalue  Vj. 

Therefore  we  can  express  the  random  vector  x  as  a  linear  combination  of 
the  eigenvectors  of  the  correlation  matrix  of  x.  Once  we  have  computed  the 
eigenvalues  and  eigenvectors,  we  will  extract  the  dominant  ones  and  use 
them  to  extend  the  original  data  sequence  by  appending  this  new  sequence  to 
the  original  sequence.  All  test  sequences  are  extrapolated  to  a  length  of  128 
points. 

B.  RESULTS 

The  first  set  of  test  data  consists  of  two  sinusoids  with  a  10  Hz  (fs  =  1000 
Hz)  difference  in  frequency  and  0  degree  relative  phase  angle.  By  examining 
the  eigenvalues  of  the  correlation  matrix,  we  observe  that  there  are  two 
dominant  eigenvectors.  This  corresponds  to  the  fact  that  we  have  two 
sinusoids  present.  Because  there  are  two  dominant  eigenvectors  present,  we 
use  the  linear  combination  of  the  two  most  dominant  eigenvectors  for  the 
extrapolation.  Therefore  in  Eq.  4.2,  N  is  set  to  2  and  the  resulting  sequence  is 
appended  to  the  original  data.  The  autocorrelation  matrix  in  Eq.  4.2  is 
computed  by  using  the  biased  autocorrelation  method.  Figure  37  shows  the 
plot  of  the  extrapolated  data  sequence  and  the  resulting  periodogram.  Two 
distinct  peaks  are  present  but  they  are  slightly  off.  The  SNR  is  15  dB. 
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Figure  37.  Eigenvector  Extrapolation  of  2  Sine  Waves,  10  Hz  Apart,  0®  Initial 
Phase  Angle,  SNR  =  15  dB,  2  Eigenvectors  Used;  and  Resulting  Periodogram 
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Figure  38  shows  the  results  when  the  initial  relative  phase  angle  between 
the  two  sine  waves  is  changed  to  45  degrees.  Once  again  two  distinct  spectral 
peaks  are  present  but  they  are  slightly  off  from  the  true  values.  Figure  39 
shows  the  results  when  the  relative  phase  angle  is  changed  to  90  degrees.  As 
in  the  two  prior  cases,  the  spectral  peaks  are  slightly  off  from  their  true 
values. 

Even  though  the  spectral  peaks  are  not  absolutely  accurate,  they  remain 
fairly  consistent  in  even  lower  SNRs.  Figure  40  shows  the  results  of  the  two 
sine  waves  (initially  in  phase,  10  Hz  apart)  in  a  SNR  of  -3  dB.  The  spectral 
peaks  are  still  well  defined.  When  the  SNR  was  lowered  to  -6  dB,  the 
periodogram  was  unable  to  distinguish  the  two  peaks.  Figure  41  shows  these 
results. 

Figure  42  shows  the  results  of  the  eigenvector  extrapolation  on  the  KM 
sequence.  N  was  chosen  to  be  3  as  we  have  three  sinusoids  present  in  this 
case.  Three  main  spectral  peaks  are  present  but  one  of  the  estimates  is  slightly 
off. 
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Figure  38.  Eigenvector  Extrapolation  of  Two  Sine  Waves,  lOHz  Apart,  45° 
Initial  Phase  Angle,  SNR  =  15  dB,  2  Eigenvectors  Used;  and  Resulting 

Periodogram 
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Figure  39.  Eigenvector  Extrapolation  of  2  Sine  Waves,  10  Hz  Apart,  90°  Initial 
Phase  Angle,  SNR  =  15  dB,  2  Eigenvectors  used;  and  Resulting  Periodogram 
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magnitude  in  dB  amplitude 


Figure  40.  Eigenvector  Extrapolation  of  2  Sine  Waves,  10  Hz  Apart,  0°  Initial 
Phase  Angle,  SNR  =  -3  dB,  and  Resulting  Feriodogram 
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Figure  41.  Eigenvector  Extrapolation  of  2  Sine  Waves,  10  Hz  Apart,  0°  Initial 
Phase  Angle,  SNR  =  -6  dB,  and  Resulting  Periodogram 
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PEIIIODOGRAM  OF  EXTRAPOLATED  SEQUENCE 


Figure  42.  Eigenvector  Extrapolation  of  KM  Data,  3  Eigenvectors  Used,  and 

Resulting  Periodogram 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


In  this  thesis,  data  extrapolation  techniques  to  enhance  the  spectral 
resolution  of  the  periodogram  are  presented.  The  three  methods  of 
extrapolation  are:  linear  predictive  filtering,  modeling  the  data  as  the 
impulse  response  of  a  filter  (Prony's  method),  and  eigenvector  extrapolation 
using  the  Karhunen-Lo4ve  expansion. 

The  most  promising  of  the  three  methods  is  the  linear  predictive  filtering 
extrapolation  using  the  prediction  coefficients  obtained  from  the  modified 
forward-backward  method.  This  technique  resulted  in  spectral  estimates 
which  were  equivalent  to  or  better  than  those  obtained  via  the  autoregressive 
method.  Frequencies  as  closely  spaced  as  2  Hz  have  been  resolved,  with  an 
appropriate  SNR,  where  normally  resolution  of  no  better  than  12.5  Hz  would 
be  obtained  using  a  periodogram  of  the  nonextrapolated  sequence.  This 
method  performed  very  well  in  low  SNRs. 

In  a  high  SNR  environment,  the  periodograms  resulting  from  Prony's 
method  of  extrapolation  gave  good  results.  Performance  was  not  as 
promising  in  lower  SNRs.  Model  order  selection  in  Prony's  method  proved 
to  be  very  important  in  determining  the  resolution  of  the  resulting 
periodogram. 

None  of  the  methods  covered  required  very  long  computer  run  times.  A 
typical  run  lasted  approximately  1  to  2  minutes  on  a  286  based  personal  com¬ 
puter,  depending  on  the  extrapolation  length.  As  discussed  earlier,  care  must 
be  taken  when  extrapolating  the  data  sequence  to  ensure  that  it  does  not  grow 
without  bound.  This  has  devastating  effects  on  the  spectrum  estimation. 
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In  this  thesis  we  have  shown  that  data  extrapolation  is  a  viable  means  for 
increasing  the  resolution  of  the  periodogram.  Further  research  should 
concentrate  on  other  data  extrapolation  techniques  for  better  and  more  robust 
spectral  resolution. 
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APPENDIX  A.  OTHER  METHODS  OF  AR  PARAMETER  ESTIMATION 


A.  AUTOCORRELATION  METHOD 

In  the  autocorrelation  method  or  the  Yule-Walker  approach,  the  AR 
parameters  are  estimated  by  minimizing  an  estimate  of  the  prediction  error 
power 


1  “ 


n=-oo 


4”]+  ^4^]x[n-k] 

ik=l 


(A.1) 


The  samples  of  x[n]  which  are  not  observed,  those  not  in  the  range  0Sn<N-l, 
are  set  equal  to  zero.  To  minimize  the  prediction  error  power,  we 
differentiate  Eq.  A.l  with  respect  to  the  real  and  imaginary  parts  of  the  a[k]s. 
This  yields 


^  n=-«>V  > 


P 

I 

n=-oov  k-\ 

In  matrix  form  this  set  of  equations  becomes 


t*[«-/]  =  0,  /  =  1,2,...P. 


(A.  2) 


'»[i) 

rxx[^]  ^xxlO]  •" 

>•«[-(? -2)] 

a[2] 

=  — 

'iipi 

/xx[P-l]  - 

rxt[0] 

.«[3]. 

— 1 

a. 

■ 

_ 1 

(A.  3) 


where 


,  N-l-ik 

-fj  +  forfc  =  0,l,...,P 

^  n=0 


(A.4) 


for  ic  =  -( p  - 1),  -(p  -  2), . . . ,  -1 


65 


which  is  recognized  as  the  biased  autocorrelation  function  (ACF)  estimator. 
The  AR  parameters  or  prediction  coefficients  are  easily  solved  for  from  Eq. 
A.3. 

* 

The  matrix  in  Eq.  A. 3  is  hermitian  (rxx[-kl  =  ^xxM)  Toeplitz,  and 

furthermore  can  be  shown  to  be  positive  definite.  As  such,  the  resulting 
estimated  poles  are  guaranteed  to  be  within  the  unit  circle.  This  property 
ensures  that  the  prediction  coefficients  calculated  from  Eq.  A.3  will  give  stable 
extrapolations.  This  is  very  important  because  if  the  extrapolation  grows 
without  bound,  it  will  have  devastating  effects  on  the  spectrum  estimation. 

The  autocorrelation  method  extrapolation  was  found  to  produce  poorer 
resolution  of  the  spectral  estimates  than  the  other  extrapolation  techniques 

I 

used.  The  length  of  the  extrapolated  sequences  in  this  section  is  320  points. 

‘  Figure  43  shows  the  resulting  periodogram  of  the  autocorrelation  method 
extrapolation  and  th'e  respective  AR  spectrum  estimation  of  the  KM  data. 
The  extrapolation  gave  a  periodogram  which  correctly  identifies  the  3 
sinewaves  present  but  the  spectral  peaks  are  not  as  well  defined  as  in  the 
periodogram  obtained  from  the  MFBLP  extrapolation  technique.  The  AR 
spectrum  estimation  resulting  from  the  autocorrelation  method  is  unable  to 
distinguish  the  two  closely  spaced  sinusoids  at  0.2  and  0.21. 

When  the  autocorrelation  method  extrapolation  was  applied  to  a  data 
sequence  consisting  of  two  sine  waves  embedded  in  white  noise  (SNR=15  dB), 
the  resulting  periodogram  was  unable  to  distinguish  the  two  closely  spaced 
sinusoids  (5  Hz  apart,  fs  =  1000  Hz).  The  AR  spectrum  estimation  also  failed. 
Figure  44  shows  the  results. 


Figvire  44.  Periodogram  Resulting  from  Extrapolation  of  2  Sine  Waves,  5  Hz 
Apart,  SNR  =  15  dB,  0°  Initial  Phase  Angle  and  Respective  AR  Spectrum 

Estimation 
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B.  COVARIANCE  METHOD 

In  the  covariance  n\ethod,  the  AR  parameters  are  estimated  by 
minimizing  the  estimate  of  the  prediction  error  power 


1  V 

^  N-P  ^ 


n=P 


'^a[k]x[n-k] 
k=\ 


(A.5) 


Note  that  the  only  difference  between  the  covariance  method  and  the 
autocorrelation  method  is  the  range  of  summation  in  the  prediction  error 
power  estimate.  The  minimization  of  Eq.  A.5  yields 


••• 

[41]] 

^xx[l'0] 

^xx[2/ll  ••• 

42] 

=  — 

^xx[2/0] 

Sxx[P>'^]  ^^xx[^/2]  — 

a[Pl 

.^xx[^'^]. 

where 

* [”  -  ;>[«  -  (A.7) 

^  ^  n=P 

The  correlation  matrix  in  Eq.  A.6  is  positive  semidefinite  and  not 
Toeplitz.  As  such,  the  estimated  poles  using  the  covariance  method  are  not 
guaranteed  to  lie  within  the  unit  circle.  Therefore  the  prediction  coefficients 
are  not  guaranteed  to  give  bounded  extrapolations  although  in  practice  the 
extrapolations  are  genercilly  bounded. 

Results  from  the  covariance  method  extrapolation  were  better  than  the 
autocorrelation  method  extrapolation.  Once  again  the  data  is  extrapolated  to 
a  length  of  320  points.  Figure  45  shows  the  periodogram  of  the  covariance 
method  extrapolation  and  the  respective  AR  spectrum  estimation  of  the  KM 
data.  The  periodogram  from  the  covariance  method  extrapolated  data  se¬ 
quence  has  well  defined  spectral  peaks  and  correctly  identifies  the  3  sinusoids 
present  The  AR  spectrum  estimation  also  identifies  the  3  sinusoids. 
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Figure  45.  Periodogram  of  Covariance  Method  Extrapolated  KM  Data  and 
Respective  AR  Spectrum  Estimate 
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When  the  covariance  method  extrapolation  was  applied  to  a  data 
sequence  consisting  of  2  sine  waves  embedded  in  white  noise  (SNR=15  dB), 
the  resulting  periodogram  was  unable  to  distinguish  the  two  closely  spaced 
sinusoids  (5  Hz  apart).  The  AR  spectrum  estimation  also  failed.  Figure  46 
shows  the  results. 

C  BURG'S  METHOD 

In  contrast  to  the  autocorrelation  and  covariance  methods  which  estimate 
the  AR  parameters  directly,  the  Burg  method  estimates  the  reflection 
coefficients  and  then  uses  the  Levinson  recursion  to  obtain  the  AR  parameter 
estimates. 

Burg's  optimization  criterion  is  to  minimize  both  the  forward  prediction 
error  ep(n)  and  the  backward  prediction  error  ep{n): 

P  =  77^  +  ep[nf  ]  (A.8) 

where 

ep[n]  =  jc[n]  +  a[l]x[n  - 1]  +  fl[2]x[n  -  2]+...+a[P]z[n  -  P] 
ep[n]  =  x[n  -  P]  +  fl[l]x[n  -  P  + 1]  +  a[2]x[n  -  P  +  2]+...+fl[P]x[n]. 

The  computional  steps  are  summarized  below  [Ref.  2:pp.  228-231] 

1.  Initialize  by  setting 

=  =  forO$n<N-l 


and  the  order  mean  square  error  as 


2.  Compute  the  reflection  coefficient: 


N-\ 

-  «=P _ 

~  N-i  ^ 

X^P-lW 

n=P 

3.  Compute  ap[k]  using  the  Levinson  recursion  given  by: 


1 

1 

0 

1 

flp[l] 

ap[2] 

«P-l[l] 

‘*P-l[2] 

-yp 

flp-i[E-i] 

ap_l[P-2] 

_ap[P] 

flp.l[P-l] 

0 

'^P-lW 

1 

4.  Compute  ep[n]  and  ep[n]  for  P  <  n  <  N-1 


(A.10) 


(A.11) 


ep[«]  =  ep-i[”]-  yp^P-i[”] 

ep[n]  =  gp_|[n-l]-  ypep_i[n].  (A.  12) 

5.  Update  the  mean-square  error  as  follows 

Ep  =(l-yp)^Ep_l. 

6.  Continue  the  iteration  until  P  equals  the  model  order. 

Figure  47  shows  the  results  using  Burg’s  method.  The  data  is  extrapolated 
to  320  points.  The  periodogram  of  the  extrapolated  KM  data  resolves  the  3 
sinusoidal  components,  however  one  of  the  spectral  components  is  slightly 
off.  In  addition  it  seems  to  be  prone  to  line  splitting.  The  respective  AR 
spectral  estimation  has  similar  results.  The  model  order  is  16. 
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Figure  47.  Periodogram  of  Extrapolated  Sequence  and  Respective  AR 
Spectrum  Estimate  of  KM  Data  Using  Burg's  Method 
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Figure  48  shows  the  PSDs  of  a  sequence  consisting  of  2  sine  waves 
embedded  in  white  noise  (SNR=15  dB)  which  are  5  Hz  apart  (fs  =1000  Hz). 
Neither  the  periodogram  of  the  extrapolated  sequence  nor  the  AR  spectrum 
estimation  could  distinguish  the  two  closely  spaced  sinusoids.  The  model 
order  once  again  is  16  and  the  data  is  extrapolated  to  320  points. 
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APPENDIX  B.  SIMULATION  RESULTS  FROM  PRONVS  METHOD 


In  Chapter  HI,  the  model  order  for  the  numerator  and  denominator  in  Eq. 
3.2  was  chosen  to  be  16.  If  too  low  an  order  is  chosen,  there  is  not  enough 
resolution  in  the  periodogram  of  the  resulting  extrapolated  sequence.  Figure 
49  shows  the  results  of  a  noncausal  extrapolation  and  the  resulting 
periodogram  for  the  KM  data  when  the  model  order  is  8.  The  spectral 
component  at  0.1  fraction  of  the  sampling  frequency  is  not  resolved  and  the 
two  closely  spaced  components  at  0.2  and  0.21  are  incorrect.  The  data 
extrapolation  length  is  noted  on  each  figure. 

If  the  order  is  too  high,  spurious  spectral  peaks  arise.  This  can  be  seen  in 
Figvu'e  50  which  is  the  noncausal  extrapolation  and  resulting  periodogram  of 
two  sine  waves  embedded  in  white  noise  (SNR=15  dB)  which  are  10  Hz  apart 
(fs  =  1000  Hz).  The  two  frequency  components  are  properly  represented  but 
an  incorrect  spectral  component  at  approximately  0.25  arises. 

When  a  causal  or  strictly  forward  extrapolation  was  tried,  the  results  were 
not  as  promising  as  the  noncausal  extrapolation.  Figure  51  shows  the  results 
of  a  causal  extrapolation  and  the  resulting  periodogram  of  the  KM  data.  The 
model  order  is  16.  The  periodogram  of  the  causal  extrapolation  was  unable  to 
resolve  the  two  closely  spaced  sinusoids  at  0.2  and  0.21.  Figure  52  shows  the 
results  of  a  causal  extrapolation  and  the  resulting  periodogram  of  two  sine 
waves  embedded  in  white  noise  (SNR=15  dB)  and  10  Hz  apart  (fs  =  1000  Hz). 
The  periodogram  from  the  causal  extrapolation  was  unable  to  resolve  the  two 
frequencies. 
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PRONYS  ESTIMATE  OF  DATA  SEQUENCE  -  ORDER  8  NON-CAUSAL 


Figure  50.  Noncausal  Extrapolation  and  Resulting  Periodogram  of  2  Sine 

Waves 
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PRONYS  ESTIMATE  OF  DATA  SEQUENCE  -  ORDER  16  CAUSAL 


Figure  51.  Causal  Extrapolation  and  Resulting  Feriodogram  of  KM  Data 
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APPENDIX  C.  RESULTS  FROM  MISCELLANEOUS  SEQUENCES 

The  four  main  types  of  extrapolation  techniques  in  this  thesis  (modified 
covariance,  modified  forward-backward,  Prony,  and  eigenvector  method) 
were  tried  on  a  noise  only  sequence.  Figures  53  through  56  show  the 
respective  results  of  the  periodogram  of  the  extrapolated  sequences  and  the 
AR  spectrum  estimations.  The  noise  sequences  are  extrapolated  to  320  points 
prior  to  the  calculation  of  the  periodogram.  Figure  57  shows  the  periodogram 
of  the  nonextrapolated  noise  sequence. 

The  MFBLP  method  and  Prony's  method  were  also  tried  on  a  test 
sequence  consisting  of  a  single  exponentially  decaying  sine  wave  with  an 
SNR  of  15  dB.  The  sine  wave  decayed  to  10%  of  its  maximum  value  within 
the  original  64  points.  Figure  58  shows  the  resulting  periodogram  from  the 
MFBLP  method  extrapolated  sequence  and  the  respective  AR  spectrum 
estimation.  Both  methods  correctly  identified  the  frequency  of  the  sine  wave. 
Figure  59  shows  the  Prony's  method  extrapolated  sequence  and  the  resulting 
periodogram.  Once  again  the  periodogram  correctly  identified  the  frequency 
present. 
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MODIFIED  COV  EXTENSION.  ORDER  16,  NOISE  SEQUENCE  ONLY 


AR  SPEa’RUM  ESTIMATION 


Figure  53.  Periodogram  of  FBLP  Extrapolation  of  Noise  and  Respective  AR 

Spectrum  Estimate 
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MODIFIED  FDLP  EXTENSION,  ORDER  48,  NOISE  SEQUENCE  ONLY 


AR  SPECTRUM  ESTIMATION 


Figure  54.  Periodogram  of  MFBLP  Extrapolation  of  Noise  and  Respective  AR 

Spectrum  Estimate 
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Figure  56.  Eigenvector  Extrapolation  of  Noise  and  Resulting  Feriodogram 
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Figure  57.  Periodogram  of  the  Nonexlrapolated  Noise  Sequence 


87 


magnitude  in  dB  amplitude 


PRONYS  ESTIMATE  OF  DATA  SEQUENCE  -  ORDER  16  NON-CAUSAL 


PRONYS  METHOD  •  ORDER  16  NON-CAUSAL 


Figure  59.  Prony's  Extrapolation  of  Exponentially  Decaying  Sine  Wave  and 
Respective  AR  Spectrum  Estimation 
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APPENDIX  D.  COMPUTER  SIMULATION  CODE 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  backfor.m 
%% 

%%  This  program  extrapolates  the  data  using  forward- 
%%  backward  linear  prediction.  The  AR  spectral 
%%  estimation  is  also  calculated. 

%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


load  kmdat.mat; 

%load  robd.mat; 
load  coeff.mat; 

a=coeff ; 

%predf=robd; 

predf=kmdat; 

freql=. 1 ; 
freq2=. 2 ; 
freq3=. 21; 

ord=16 ; 
ext*256; 
pad=1024 ; 

%  Forward  prediction 


%  load  data  sequence 

%  load  prediction/AR 
%  parameters 


%  model  order 
%  amount  extrapolated 


n=65; 

for  i=l:ext; 

p(i)*linpred(predf ,n,a,ord) ; 

predf (n)=p(i) ; 

n=n+l; 

end 

%  Backward  prediction 

flip=flipud(kmdat) ;  %  extrapolate  backwards 

m=65; 

for  i=l:ext; 

p(i)=linpred(flip,m,a,ord) ; 
flip{m)*p(i) ; 
m®m+l ; 

end 

predb=flipud(flip) ; 
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total=[predb(l:ext) ;predf ] ; 
clear  predf  predb  flip  p  a  kmdat 


%plot (total) , title ( 'EXTRAPOLATED  SEQUENCE') 

%xlabel ( ' sample ' ) ,  ylabel ( ' amplitude ' ) 

%meta  bflp 
%r'use 

k=length (total) 
window=hamraing(k) ; 
f inham=window  .*  total; 
clear  window  total 
pow=fft(f inham, pad) ; 
n=length(pow) ; 
pyy=pow  .*  conj  (pow) 

f=0.5*(0;n/2)/(n/2)  ; 

pyy(n/2+2:n)=[]; 

pownorm= ( 1/max (pyy) )  .*  pyy; 

powlog=10  .*  logic (pownorm) ; 

clear  pow  pyy  pownorm  k 

fl=[freql  -100;  freql  100]; 
f2=[freq2  -100;  freq2  100]; 
f3=[freq3  -100;  freq3  100]; 
axis((0  .5  -60  0])  ; 

a=[l;coeff ] ; 

sar=abs(fft (a,pad) )  .'2; 
sar=  1.0  ./  sar; 
sarnorm= ( 1/max ( sar) )  .*  sar; 
sarlog=  10  .*  loglO (sarnorm) ; 
sarlog=sarlog(l;n/2  +  1) ; 

plot(f,powlog,'-',fl(:,l),fl(;,2),'  — •,f2(:,l),f2(:,2)) 

title ( 'PERIODOGRAM  OF  EXTRAPOLATED  SEQUENCE') 

ylabel ( 'RELATIVE  PSD  (dB)') 

xlabel (' FRACTION  OF  SAMPLING  FREQUENCY') 

meta  bflp3 

pause 

plot(f,sarlog,  • ,  f  1  ( : ,  1)  ,  f  1  ( : ,  2) ,  •  — ' ,  f2  ( : ,  1)  ,  f2  ( : ,  2) ) 

title ('AR  SPECTRUM  ESTIMATION') 

xlabel ( 'FRACTION  OF  SAMPLING  FREQUENCY') 

ylabel ( 'RELATIVE  PSD  (dB) ' ) 

meta 

pause 

axis; 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function  lp=linpred(x,n,a,ord) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  linpred  computes  the  nth  element  of  a  n-1  sequence  using 
%  simple  linear  prediction,  x-data  sequence;  n=element  to  be 
%  predicted;  a=filter  coefficients;  ord=order  of  predictor. 

i=l; 

for  k=l;ord; 

y (i)=a(k) *x(n-k) ; 
i=i+l ; 

end 

lp=-sum(y) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%  mfblp.m 
%% 

%%  This  program  computes  the  prediction/AR  parameters 
%%  using  the  modified  forward/backward  method. 

%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%load  kmdat.mat; 
load  robd.mat; 

%x=kmdat ; 
x=robd ; 

n=length(x) ;  %  number  of  data  points 

m=96; 

np=n-m ; 

for  i=l:np;  %  load  data  matrix  A 

for  k=l:m; 

A(i,k)=x(i+ra-k) ; 

A(i+np,k)=x(i+k) ; 

end 

end 

%[u,s,v]»svd(A) ;  %  singular  value  decomposition 

[ eigvec , eigval ] =eig (A ' *A) ; 


k=4; 


%  number  of  sinusoids  present 
%  in  signal,  assumed  known 
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%  eigenvalues  of  A' A 


eigval=diag(eigval) ; 

%eigval=diag(s)  .“2; 
w=zeros (in,  1)  ; 

for  i=l:n-in; 

b(i,l)=x(in+i)  ;  %  2(n-Tn)-by-l  desired 

b(i+n-in,  l)=x(i)  ;  %  response  vector 

end 

for  i=l:k; 

w=w  +  (eigvec( ;  ,in-i+l)  ./  eigval  (m-i+l)  *eigvec( :  ,in-i+l) 

end 


%for  i=l:k; 

%  w=w  +  (v(:,i)  ./  eigval (i)  *  v(:,l)'  *  A'  *  b) ; 
lend 

1  erase  coeff.mat 
diary  coeff.mat 
w=-w 

diary  off 
!q  coeff.mat 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%  pron.m 
%% 

%%  This  program  finds  an  HR  filter  with  a  prescribed 
%%  time  domain  response  using  Prony's  Method  for  filter 
%%  design.  The  impulse  response  of  that  filter  is  then 
%%  computed  and  a  periodogram  is  computed  to  estimate 
%%  the  spectrum  of  the  data.  The  data  is  extrapolated 
%%  by  increasing  the  impulse  duration. 

%% 

%%%%%%%%%%%%l%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


load  kmdat.mat; 

h=kmdat;  %  time  domain  impulse  response 

hf=flipud(h) ; 


nb=16; 
na=16; 
zpad=1024 ; 


%  order  of  numerator 
%  order  of  denominator 


[b,a]=prony (h,nb,na)  ;  %  HR  filter  coefficients 

[d,c]=^ony  (hf  ,nb,na)  ; 

n=128; 

x=tl  zeros(l,n-l) ] ; 
y=filter(b,a,x) ;  ynew=y (65:n)  ' ; 

yf=filter (d,c,x) ;  ynewf=flipud(yf (65:n)  ')  ; 
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